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1.  Introduction 

The  study  of  nonlinear  acoustic  phenomena  in  dissipative 
fluids  has  been,  and  remains,  an  important  branch  of  fluid 
mechanics.  In  the  1950s,  Mendousse  [1],  and  later  Lighthill  [2], 
showed  that  under  certain  assumptions,  the  essential  one  being 
that  fluctuations  in  the  field  variables  are  sufficiently  small  in 
amplitude,  Burgers’  equation  provides  an  approximate  description 
of  “finite  amplitude”  acoustic  waves  in  dissipative  perfect  gases; 
see  also  Ref.  [3]  and  those  therein.  Since  that  time  researchers 
have,  in  an  effort  to  obtain  ever  more  accurate,  but  still  tractable, 
descriptions  of  acoustic  propagation  in  gases  and  liquids,  derived 
a  number  of  other  so-called  weakly-nonlinear  model  equations, 
i.e.,  nonlinear  approximations  to  the  compressible  Navier-Stokes 
equations. 
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The  primary  aim  of  the  present  study  is  to  determine,  for 
physically  realistic  parameter  regimes,  which  of  the  following 
bi-directional,  weakly-nonlinear  models  best  approximates  the 
compressible  ID  Navier-Stokes  equations  under  the  travelling 
wave  assumption. 

Lighthill-Westervelt  (LW)  equation  [4,5]: 

cjcpxx  -  4>tt  +  8c~24>m  =  pc~2dt{4>t)2.  (1) 

Blackstock-Lesser-Seebass-Crighton  (BLSC)  equation  [6,7]: 

[cl  -  208  -  l)0t]0xx  -  <Ptt  +  S(ptxx  =  dt(4>x)2-  (2) 

Kuznetsov’s  equation  [8,9]: 

cl<pxx  -  cptt  +  8<t>txx  =  dtlicpx)2  +  c~2(P  -  1  )(<Pt)2].  (3) 

Rasmussen-Sorensen-Gaididei-Christiansen  (RSGC)  equation 

[10,11]: 

(cl  -  4>t)<t>xx  -  4>tt  +  =  3 1  |W  +  C~2  (p  -  0  (< Pt )2  . 

(4) 

These  four  models  are  based  on  the  assumption  that  the  propa¬ 
gation  medium  is  a  classical  thermoviscous  fluid5  whose  equilib¬ 
rium  (or  undisturbed)  state  is  characterised  by  q  =  Qe,  p  = 


5  By  which  we  mean  a  Newtonian  fluid  whose  transport  coefficients  are 
constants,  and  within  which  the  flow  of  heat  is  governed  by  Fourier’s  law;  see, 
e.g.,  [12,  Chap.  10]. 
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pe,  u  =  (0,  0,  0),  0  =  0e ,  and  rj  =  rje,  where  q,  p,  u,  0  and  r] 
denote  the  mass  density,  thermodynamic  pressure,  velocity,  ab¬ 
solute  temperature,  and  specific  entropy,  respectively.  Addition¬ 
ally,  all  material  parameters,  as  well  as  all  quantities  with  an  “e” 
subscript,  are  regarded  as  constant;  the  velocity  vector  is  taken  as 
u  =  (u(x,  t),  0,  0),  i.e.,  planar  flow  along  the  x-axis  is  assumed; 
u  =  0X,  where  0  =  0(x,  t)  denotes  the  velocity  potentially  virtue 
of  the  fact  that  the  irrotationality  condition  V  x  u  =  0  is  identi¬ 
cally  satisfied  here;  and  all  body  forces  are  neglected.  Moreover, 
ce(>0),  the  adiabatic  sound  speed,  denotes  the  sound  speed  in  the 
undisturbed  fluid;  is  known  as  the  coefficient  of  nonlinearity  [13], 
where  e  (1,7)  for  nearly  all  fluids  under  ordinary  conditions; 
and 

5  =  v[4/3  +  ^B/^  +  (y-l)/Pr],  (5) 

termed  the  diffusivity  of  sound  by  Lighthill  [2],  is  a  positive 
constant  [14,  Sect.  4.13].  In  Eq.  (5),  the  constants  /x(> 0)  and  /xB 
respectively  denote  the  shear  and  bulk  viscosities,  where  we  note 
that  /xB  is  generally  nonzero  [14,  p.  21];  v  =  /x/ge  is  the  kinematic 
viscosity;  the  adiabatic  index  is  defined  as  y  :=  cp/cv ,  where  the 
constants  cp  >  cv  >  0  respectively  denote  the  specific  heats  at 
constant  pressure  and  volume  [  14, 15  ] ;  and  Pr  =  v/k  is  the  Prandtl 
number,  where  k  is  the  thermal  diffusivity. 

The  approach  we  adopt  here  is  to  first  obtain,  wherever  possi¬ 
ble,  analytical  results,  i.e.,  exact  and/or  approximate  solutions,  and 
then  turn  to  numerical  methods  for  the  purposes  of  comparison 
and  illustration.  The  present  study  could  be  regarded  as  a  follow- 
on  to  those  presented  in  Refs.  [11,  Sect.  5.2]  and  [16],  the  aims  of 
which  were  to  compare  several  lossless  weakly-nonlinear  models 
with  the  Euler  equations. 

In  drawing  this  section  to  a  close,  we  feel  it  necessary  to 
point  out  the  following.  The  four  PDEs  we  consider  are  all,  as  a 
survey  of  the  acoustics  literature  reveals,  well-established,  weakly- 
nonlinear  approximations  of  the  ID  Navier-Stokes  equations. 
They  have  been  derived,  as  alluded  to  earlier,  primarily  so  that 
researchers  might  have  simplified  (i.e.,  analytically  tractable) 
alternatives  to  the  latter;  in  other  words,  as  single-equation 
models  that  are,  nonetheless,  capable  of  capturing  the  nonlinear 
phenomena  exhibited  by  Newtonian  fluids  undergoing  irrotational 
compressible  flow.  To  be  clear,  our  purpose  here  is  not  to 
attempt  to  improve  upon  Eqs.  (l)-(4)  by  including  higher 
order  nonlinearities;  indeed,  doing  so,  assuming  all  requisite 
higher-order  terms  in  the  constitutive  relations  are  known 
unambiguously  [15,  Sect.  1.1.5],  is  not  only  counter  to  the  reason(s) 
these  PDEs  were  originally  derived,  but  negates  their  usefulness  as 
models  by  making  them  more  complicated,  possibly  to  the  point 
where  the  analytical  challenges  they  would  then  pose  rivaled  those 
of  the  (ID)  Navier-Stokes  equations  themselves. 

Remark  1.  In  the  case  of  perfect  gases,  ce  =  V YPe/Qe  and  can 
be  expressed  in  terms  of  y,  specifically,  by  /3  =  (y  +  l)/2. 

2.  Governing  equations 

All  four  weakly  nonlinear  models  are  based  on  the  continuity 
and  momentum  equations,  which  in  ID  reduce  to 


Qt  +  UQx  +  QUx  =  0,  (6) 

g(ut  +  uux)  =  —px  +  Uxx »  (7) 

respectively,  recalling  that  all  body  forces  have  been  neglected,  as 
well  as  the  linearised  ID  energy  equation  [17,  p.  21] 

Qe^eht  =  K0XX  =  XK  Qe0e(ptxx,  (8) 


and  the  assumption  that  the  equation  of  state  can  be  taken 
as  [15,  pp.  581-583] 


P  =  Pe  +  QeCe 


s  +  (/3  —  l)s  + 


y-i 

xcl 


(r]  -  Tie) 


(-1  «  s  «  1),  (9) 

i.e.,  as  a  truncated  Taylor  expansion  of  the  general  constitutive 
equation  p  =  p(q,tj).  Here,  s  =  (q  —  Qe)/Qe  denotes  the 
condensation  and  the  positive  constants  I<  and  x  represent  the 
thermal  conductivity  and  the  thermal  coefficient  of  volume 
expansion,  respectively,  where  it  should  be  noted  that  I<  =  KCpge. 


Remark  2.  As  alluded  to  earlier,  since  they  are  based  on  Eq.  (9), 
Eqs.  (1)— (4)  describe  acoustic  waves  in  both  liquids  and  gases, 
provided  of  course  that  fluctuations  in  q,  u,  etc.,  about  their 
equilibrium  state  values  are  sufficiently  small. 


3.  TWSs:  weakly-nonlinear  models 

In  this  section  we  review  and  examine  exact  TWSs  of  the 
four  weakly-nonlinear  models.  To  this  end,  we  begin  by  recasting 
Eqs.  (l)-(4)  in  terms  of  the  following  dimensionless  quantities: 

0°  =  0/(t/L),  u°  =  u/U,  x°  =  x/L,  t°  =  t  (ce  /  L) .  (10) 

Here,  the  positive  constants  U  and  L  denote  a  characteristic  speed 
and  length,  respectively,  and  for  typographical  convenience  all 
diamond  superscripts  will  henceforth  be  omitted  but  are  to  remain 
understood.  For  future  reference  we  note  that  e  =  U/ce  is  the 
Mach  number,  where  e  1  is  henceforth  assumed;  Red  =  ceL/8 
denotes  a  Reynolds  number;  §(■)  :=  x  —  (■)£  is  the  wave  variable; 
a  prime  denotes  d/d§ ;  and  n  is  an  index  which  assumes  the  values 
n  =  0,  1,2,3. 

And  lastly,  since  Eqs.  (l)-(4)  are  all  invariant  under  the 
transformation  x  i->  — x,  we  henceforth  assume,  without  loss  of 
generality,  only  right-running  waves. 


3.1.  The  LW  and  BLSC  equations 


In  dimensionless  form,  Eqs.  (1)  and  (2)  become,  respectively, 

0xx  -  <ptt  +  (RedrVttt  =  epM4>t)2.  (11) 

[1  -  2e(P  -  1  )</>(]</>„  -  4>tt  +  (RedrVixx  =  2 e<t>x<t>xt.  (12) 

Introducing  the  ansatzes  u(x,  t)  =  F($(vn))  and  0(x,  t)  = 
P($(vn)),  where  F  :=  the  constants  un(>0)  denote  the  speeds 
of  the  travelling  waves,  and  in  this  subsection 


JO,  LW, 
[1,  BLSC, 


(13) 


integrating  each  once  w.r.t  §,  and  then  imposing  the  asymptotic 
condition6  F  — ►  0  as  §  — ►  oo,  Eqs.  (11)  and  (12)  reduce  to 


(Red)-1!7'  -  (Lpf)  F  =  e/3F2  (n  =  0,  1),  (14) 

i.e.,  the  associated  ODEs  of  the  former  and  latter  are  Bernoulli 
equations. 

Integrating  these  ODEs  and  then  imposing  the  usual  wavefront 
condition  /(0)  =  1/2,  the  following  Taylor  shock  solutions  are 
obtained: 


F(Hvn))  =  ^{1  -  tanh[2£ (vn)/£n]}  (n  =  0,  1).  (15) 


6 


Physically,  this  implies  that  the  fluid  at  £  =  +oo  is  in  its  equilibrium  state. 
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Here,  the  shock  thicknesses,  ln,  and  the  vn  are  given  by 


Au  — 


4 

c^RecT 


Vn 


1  H-  2  COS  [  ^  (0  T-  47T  )  J 

JTp 


n  =  0, 

(16) 

n  =  1, 


where  cos(0)  =  (2  —  2Te2p2)/2  follows  from  Cardan’s  formula 
[18].  Moreover,  the  vn  correspond  to  the  smallest,  or  only,7  positive 
root  of  the  equations 


epv3n~2n  -  v2n  +  1  =  0,  (17) 

which  arise  from  the  corresponding  ODEs  in  Eq.  (14)  as  a 
consequence  of  enforcing  the  (second)  asymptotic  condition  F  -> 
1  as  §  — >►  —  oo. 

It  is  an  interesting  fact  that,  like  those  of  the  n  =  2,  3  cases  (see 
Section  3.2),  the  wave  speed  parameter  of  the  n  =  0  case  suffers 
a  bifurcation;  specifically,  there  exists  the  critical  Mach  number 
value 

4(P)  =  P*o/P>  where  p*0  :=  «  0.385, 

such  that  v0  e  (1,^3]  fore  e  (0,  cj],  but  where  v0  e  €  for 
e  >  €q.  Here,  it  is  understood  that  taking  e  =  6q  is  possible  only 
for  p  Pq  ,  since  e  1  is  required. 

It  should  also  be  noted  that  when  c  e  (0,  ej),  the  n  =  0  case  of 
Eq.  (17)  admits  a  second  positive  root,  denoted  here  by  V0,  where 
V0  e  (vo,  oo);  however,  this  rootwas  rejected  due  to  its  unstable 
nature,  specifically,  the  fact  that  V0  ->  oo  as  e  ->  0. 

Finally,  we  observe  that  iq  >  1  for  common  fluids  under  ordi¬ 
nary  conditions,  which  follows  immediately  from  the  fact  that 


«i  =  l  +  \ep  +  \e2p2  +  ■■■  (e  «  1),  (18) 

Z  O 

and,  in  turn,  that  for  a  given  Mach  number,  v0  >  v\  on  1  <  p  < 
Pq/€,  the  interval  over  which  v0  e  R. 


Remark  3.  When  the  Mach  number  is  expressed  in  terms  of  its 
critical  value  via  the  relation  e  =  A0Cq,  where  A,0  e  (0,  1],  v0 
becomes  independent  of  p. 


Remark  4.  The  BLSC  equation  is  easily  derived  from  what  has 
come  to  be  known  as  the  Soderholm  equation  by  simply  omitting 
terms  of  (9(c2)  in  the  latter;  see,  e.g.,  Refs.  [11,19]  and  those 
therein. 


3.2.  The  Kuznetsov  and  RSGC  equations 

In  dimensionless  form,  Eqs.  (3)  and  (4),  reduce  to  the  indicated 
special  cases  of  the  following  “dummy”  PDE: 

0xx  —  0tt  (R^d)  0txx  —  c(/l  —  2 )0t0xx 

+  +  ,(19) 

where 

{2,  Kuznetsov, 

3,  RSGC, 

and  we  stress  that  this  PDE  does  not  have  a  (known)  physical 
meaning  for  values  of  n  other  than  two  and  three. 

Referring  the  reader  to  Refs.  [9, 1 1  ]  for  specifics  on  the  cases  n  = 
2,  3,  respectively,  we  observe  that  under  the  same  ansatzes  and 
conditions  introduced  in  the  previous  subsection,  the  associated 


ODEs  of  Eq.  (19),  like  those  of  Eqs.  (11)  and  (12),  are  of  the  Bernoulli 
type,  namely, 


(Redr’F'  ■ 


1  -  vt 


+  \  P  —  -n  \  Vn 


(20) 


which  when  integrated  yield  the  Taylor  shock  profiles 


F(Hvn))  =  ^{1  -  tanh[2£ (vn)/in]}  (n  =  2,  3).  (21) 

Here,  v2,3  denotes  the  smallest,  or  only,  positive  root  of  the 
corresponding  equations  nn(vn)  =  0,  where  we  have  defined 


nn(vn)  := 


P  -  nn)vn~vn  +  2€nVn  +  1 


in  =  2,  3);  (22) 


the  shock  thicknesses  are  given  by 
ll'3  Red  [e  (p  -  |n)  vl  +  |en] 

4l>2.3 


Red (^2,3  “  !) 


(n  =  2,3), 


(23) 


where  we  observe  that  u2,3  >  1;  and  the  restriction  p  /  3/2  is 
required  for  the  case  n  =  3. 

For  p  g  (1,  3/2),  a  range  which  includes  the  ^-values  of  com¬ 
mon  (e.g.,  diatomic  and  monatomic)  gases  under  ordinary  condi¬ 
tions,  it  is  not  difficult  to  establish  that  v3  e  (1,2)  and  that  this, 
the  only  positive  root  of  773  =  0  when  p  e  (1,  3/2),  is  of  multi¬ 
plicity  one. 

If,  on  the  other  hand,  we  assume  p  ^  /3|s,  where 


32 

—  «  1.185 
27 


and 

and  set  c  =  €*(P),  where 


27+7V2T 


36 


1.641, 


<(P) 


s/]3(9/3  -  8)3/2  -  (27 P2  -  36/0  +  8) 


808  -  1) 


n  =  2, 


9jP(P  -  4/3)3  -  (9 p2  -  18 p  +  6) 


(24) 


90 8  -  3/2) 

then  it  can  be  shown  that 
max(i/n) 

1  +  Vl-  303  -l)[e2*Q3)F 
3 (P  ~  l)e2*0 6) 

2  +  ^4-  18(/0  -  3/2)[€*(/0)]2 
6(P-3/2)e*(fi) 


n  =  3, 


n  =  2 


n  =  3 


(P  »  K), 


(25) 


each  of  which  is  a  root  of  multiplicity  two.  Here,  we  have  imposed 
large-/!  requirements  to  ensure  e\  3  4C  1,  as  the  weakly-nonlinear 
approximation  demands,  just  as  we  did  in  Section  3.1  to  ensure 

*o  «  1- 

When  p  >  1,3/2  and  c  e  (0,  e£3),  the  equations  772)3  =  0 
admit  two  distinct  positive  roots,  the  larger  of  which  in  each 
case,  denoted  here  as  V2;3,  we  reject  since  they  blow-up  as 
€  — >  0.  Moreover,  c|3  are  Mach  number  values  at  which  the 
corresponding  wave  speeds  exhibit  a  bifurcation:  see  Ref.  [9,  Fig.  1  ]8 
for  the  case  n  =  2.  It  is  also  noteworthy  that 

V3  <  max(i/2)  <  max(i/3)  (P  >  3/2), 

where  max(u2;3)  ->  ^3  as  p  ->►  oo;  i.e.,  max(u2)3)  ->  max(i/0)  as 
P  — >►  oo. 


7  In  this  article,  “only  positive  root”  also  refers  to  a  single  positive  root  of 
multiplicity  two. 


8  The  quantities  Vo,i,2,ec,  and  /3i  that  appear  in  Ref.  [9]  correspond  to 
max(u2),  v2,  V2,  and  respectively,  of  the  present  article. 
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Remark  5.  On  replacing  the  “small”  terms  (j)tXx  and  (0X)2  in  the 
n  =  2  case  of  Eq.  (19)  with  and  (0t)2»  respectively,  based  on 
the  (9(1)  approximations  0XX  ~  0tt  and  0X  ~  — 0t  [7],  one  obtains 
Eq.  (11),  i.e.,  the  LW  equation. 

Remark  6.  Unlike  Kuznetsov’s  equation,  the  RSGC  equation  ad¬ 
mits  a  Hamiltonian  structure  in  the  lossless  (i.e.,  Red  ->  oo) 
limit,  as  do  the  Euler  equations  to  which  this  limiting  case  corre¬ 
sponds;  see  Ref.  [11,  Sects.  2.3,  2.4],  wherein  it  is  also  shown  how 
the  BLSC  and  Kuznetsov  equations  are  derivable  from  the  RSGC 
equation. 

3.3.  Energy-rate  results 


a  dE/dt 


Seeking  further  insight  into  the  physics  captured  by  these 
weakly-nonlinear  models,  we  now  perform  an  energy  analysis  [20] 
based  on  their  TWSs.  To  this  end,  we  return  to  Eqs.  (11),  (12)  and 
(19),  multiply  each  by  0t,  and  then  integrate  over  the  real  line  w.r.t 
x.  Thereupon,  after  integrating  by  parts,  simplifying,  and  grouping 
terms,  we  are  led  to  consider 


dE 

dt 


2  -  d  [ +°° (A.  \3 
J  — ( 


-O  -  -e/3  — 
3  at 


H)£ 


(fr)3  dx, 

+oo 


-0  -  2e  ^  -  - 

/+oo 

(0t)20xx  dX 

-oo 

2  d  r+°° 


HHL 


(fry  dx, 

+oo 


n  =  0, 


(fr)  frx  dx,  n  =  1, 


(26) 


n  =  2, 


&  ~  -e  (  P  -  -  )  —  /  (fr)3  dx,  n  =  3, 


which  we  should  point  out  are  all  exact  relations  in  the 
sense  that,  apart  from  those  used  in  the  original  derivations  of 
Eqs.  (l)-(4),  no  other  approximations  have  been  employed.  Here, 
limx^±00  0X,  limx^±00  <pt ,  and  limx^±00  0tx  were  evaluated  using 
the  fact  that  the  Taylor  shocks  given  in  Sections  3.1  and  3.2  imply 
0X  =  f(x  -  vnt)  and  0t  =  -vnF(x  -  vnt),  for  n  =  0,1,  2,  3;  the 
energy  functional,  E ,  is  defined  as 


\  p+oo 

E:=Er+-  [(fr)2  +  (fr)2]  dx 

^  J —oo 


\v0,  n  =  0, 

|u„(l  +  evn),  n  =  1,2,  3, 


(27) 


where  ER  is  a  constant;  and  0  denotes  the  usual  dissipation 
integral,  namely, 


0  := 


1 

Red 


1 0t  I  (pttt  dx, 
(frt)2dx. 


n  =  0, 
n  =  1,2,3. 


(28) 


Because  their  integrands  involve  (at  most)  products  of  hyper¬ 
bolic  functions,  the  integrals  on  the  RHS  of  Eq.  (26)  can  be  evalu¬ 
ated  in  closed  form.  Omitting  the  details,  it  is  readily  established, 
using,  e.g.,  one  of  the  many  commercially  available  symbolic  soft¬ 
ware  packages,  that 


K, 

n  =  0, 

(P  -  2)vj, 

n  =  1, 

[(P  -  \)v\  -  Uv22, 

n  =  2, 

(29) 

[(P  ~  3/2)i>3  -  t 

n  =  3, 

b  dE/dt 


Fig.  1.  (a,b)  dE/dt  vs.  in  the  case  of  gases  and  liquids,  respectively,  for  e  =  0.1. 
Thick-solid:  LW (n  =  0).  Thin-solid:  BLSC(n  =  1).  Thin-broken:  Kuznetsov  (n  =  2). 
Thick-broken:  RSGC  (n  =  3). 


which  we  observe  is  independent  of  the  Reynolds  number.9  We 
also  observe  that  for  n  =  1,2,  3,  dE/dt  g  0  in  the  case  of  gases 
(i.e.,  1  <  /3  <  1.35)  and  liquids  (i.e.,  ^  >  2),  respectively,  while 
dE/dt  >  0  for  both  gases  and  liquids  when  n  =  0. 

To  illustrate  these  results,  we  have,  in  Fig.  1,  plotted  the 
expressions  for  dE/dt  given  in  Eq.  (29)  as  functions  of  p,  but  with 
the  Mach  number  fixed  at  6  =  0.1.  Therein,  we  observe  that,  for 
both  gases  and  liquids,  the  curves  for  the  n  =  1,2,3  cases  are  not 
only  similar  in  appearance  and  slope  but  are  also  grouped  rather 
closely  together,  with  all  appearing  below  (resp.  above)  the  /?-axis 
in  the  case  of  gases  (resp.  liquids).  In  contrast,  the  n  =  0  curve, 
which  is  defined  (i.e.,  real-valued)  here  in  the  case  of  liquids  only 
for  2  <  p  <  e~xpl  ~  3.849,  clearly  is  removed  from  the  others, 
always  lying  above  the  ^-axis  and  always  exhibiting  a  positive 
slope. 


Remark  7.  In  obtaining  the  expressions  given  in  Eq.  (29),  we  found 
that,  on  evaluation,  the  dissipation  integral  could  be  expressed  as 


ePv 


4—2  n 
n 


Vn(vt  -  1), 


n  =  0, 1, 
n  =  2,  3, 


(30) 


where,  as  expected  based  on  Eq.  (29),  0  in  each  case  turns  out  to 
be  independent  of  Red.  It  is  also  noteworthy  that,  under  each  of  the 
weakly-nonlinear  models,  0  is  both  strictly  positive,  as  the  physics 
of  thermoviscous  flow  demands,  and  an  increasing  polynomial 
function  of  the  vn. 


^  Meaning  that  dE/dt  is  also  independent  of  5;  see  Ref.  [12,  p.  616],  wherein 
mention  of  this  point  is  made  in  relation  to  Burgers’  equation. 
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4.  TWS:  ID  Navier-Stokes  equations 

Our  aim  in  this  section  is  to  determine,  by  analytical  methods, 
a  standard  against  which  the  TWSs  presented  in  Section  3 
can  be  meaningfully  compared/contrasted;  specifically,  to  solve 
Eqs.  (6)-(9)  under  the  traveling  wave  assumption.  While  this  task 
might  appear  all  but  impossible,  given  the  number  of  equations 
involved  and  their  mathematical  complexity,  a  kink-type  TWS 
of  this  system  can,  nevertheless,  be  obtained  by  regarding  the 
flow  as  homentropic  [14,21],  which  refers  to  the  q  =  qe  special 
case  of  isentropic  flow,  as  demonstrated  recently  by  Rasmussen 
[11,  Eq.  (3.41)]. 


41.  Ansatzes,  associated  ODE,  and  stability  results 


We  now  proceed  under  the  homentropic  flow  approxima¬ 
tion,10  which  translates  here  into  letting  K  —>  0  and  then  in¬ 
tegrating  equation  (8)  subject  to  q(x,  0)  =  rje.  To  this  end,  we 
return  to  Eqs.  (6),  (7)  and  (9)  and  use  the  third,  which  now  as¬ 
sumes  barotropic  [14]  form,  to  eliminate  p  from  the  second.  After 
recasting  the  resulting  equations  in  dimensionless  form  and  set¬ 
ting  vL  :=  v(4/3  +  /xB//x),  we  end  up  with  the  coupled  pair 

St  +  €USX  +  €UX(  1  +  s)  =  0,  (31) 

e(l  +  s)(ut  +  euux)  =  -[s  +  (fi  -  l)s2]x  +  e(R e)_1um,  (32) 

where  Re  is  a  Reynolds  number  based  on  vL.  Here,  in  the  spirit 
of  the  terminology  of  Hayes  [22,  p.  38],  we  call  14  e  [v,  8)  the 
kinematic  longitudinal  coefficient  of  viscosity;  it  represents  the 
purely  viscous  contribution  to  the  diffusivity  of  sound,  where  we 
observe  that  8  ->  vL  (from  above)  as  I<  ->  0. 

Again  seeking  right-running  waves,  we  set  u(x,  t)  =  f($)  and 
s(x,  t)  =  g($),  where  $(c)  :=  x  —  ct  and  c  is  a  positive  constant, 
and  then  substitute  these  ansatzes  into  Eqs.  (31)  and  (32).  On 
integrating  the  former  with  respect  to  and  then  using  the  fact 
that  u  =  s  =  0  in  the  equilibrium  state  to  solve  for  the  resulting 
constant  of  integration,  our  system  of  PDEs  is  reduced  to 

g  =  ef(c-ef r1  (f  //s),  (33) 

-6(1  +g)f'(c  -  ef)  =  -{g  +  (P~  1  )g2]'  +  e(R e)-1/",  (34) 

where /s  :=  cc-1  and  we  recall  that  a  prime  denotes  d/df . 

Next,  eliminating  g  between  the  last  two  equations,  integrating 
the  result  once,  and  then  enforcing  the  first  of  our  (two)  asymptotic 
conditions,  namely,  /  ->  0  as  §  ->  00,  the  associated  ODE 
for  /  turns  out  to  be  the  following  special  case  of  Darboux’s 
equation  [23]: 


(Re)-1/' 


e(fi  ~  l)/2  +/(c  -  ef)  -  c/(c  -  ef)2 
(c  -  ef)2 


(35) 


where  the  ensuing  (i.e.,  second)  constant  of  integration  is 
necessarily  zero.  Here,  we  observe  that  this,  rather  unusual  ODE 
admits  the  three  equilibrium  solutions 


f2  =  0, 


/l,3  — 


P  +  2(c2  -  1)  T  VQM 

2  6C 


(36) 


where  we  have  set  Q(c)  :=  p2  +  4(/3  —  1) (c2  —  1)  for  convenience, 
and  it  should  be  noted  that  Q(c)  is  strictly  positive.  A  phase-plane 
analysis  reveals  that  /  =  /1j2,3  are  unstable,  stable,  and  stable, 
respectively,  for  all  physically  allowable  values  of  e  and  p,  where 
we  observe  that  f2  </i  </s  <fs- 


10  Of  course,  only  lossless  fluids,  i.e.,  fluids  which  are  neither  viscous  nor  thermally 
conducting,  are  capable  of  exhibiting  homentropic  flow;  however,  as  we  will  see, 
the  impact  on  the  weakly-nonlinear  models  of  assuming  r\  =  is  simply  a  decrease 
in  the  value  of  8. 


Remark  8.  From  Eq.  (5)  it  is  clear  that  the  approximation  8  ~  14 
applies  only  to  fluids  for  which  Pr  ]£>  y  —  1  is  satisfied,  examples 
of  which  include  glycerin  [14,  Table  4.6],  pure  water  [15,  p.  583], 
and  oils  [24,  p.  80]. 


4.2.  Derivation  and  analysis  of  wave  speed  expressions 


Enforcing  now  the  second  asymptotic  condition,  i.e.,/  ->  1  as 
§  ->  —00,  we  find  that  the  speed  of  our  travelling  wave  profile 
corresponds  to  the  positive  root  of  A(c)  =  0,  where  we  have 
defined 


A{c)  :=  c3  -  2ec2  -  (1  -  e2)c  -  e(P  -  2),  (37) 

that  renders  f\  =  1.  With  regards  to  the  zeros  of  Eq.  (37)  we 
note  the  following:  (a)  for  p  >2  (i.e.,  liquids),  A(c)  =  0  has 
exactly  one  positive  root,  the  value  of  which  is  always  greater 
than  one;  (b)  for  p  <  2  (i.e.,  gases),  this  cubic  equation  has  two 
unequal  positive  roots,  at  least  one  of  which  lies  on  the  interval 
(0,  1);  and  (c)  because  it  only  arises  in  the  context  of  multiphase 
flows  (specifically,  those  involving  certain  dusty  gases  [  15,  p.  58 1  ]), 
which  of  course  are  not  described  by  the  Navier-Stokes  equations, 
the  case  p  =  2,  which  would  result  in  a  single  positive  root  of 
multiplicity  two,  will  not  be  considered  here. 

Since  we  are  only  interested  in  kink- type  solutions  of  Eq.  (35), 
we  take  c  =  cm,  where  f\  =  1  only  when  c  =  cm,  and,  of  course, 
require /(0)  e  (0,  1),  i.e., /(0)  e  (f2,/i).  Here,  using  Cardan’s 
formula  [18]  once  again,  cm,  the  largest  (or  only )  positive  root  of 
A(c)  =  0,  is  found  to  be 


Cm 


^  e  +  cos  (4)  4  +  e2  , 

2e  3  +  e2 

3  9[(27eOS-4/3)-2e3)/54  +  v/Zj1^3 

+  [(27c 08  -  4/3)  -  2e3)/54  +  s/~A  ]  V3  , 


A  <  0; 


A  >  0; 


where  A  =  A(e,/3),  the  discriminant  of  Eq.  (37),  is 


(38) 


A(e,/3) 


-4(y3  -  l)e4  +  (27 p2  -  72/3  +  44)e2  -  4 
108  ' 


and  the  angle  &  is  defined  by  the  relation 


(39) 


cos(tf)  =  €- 


21  (P  -  4/3)  -  2e2 

V(3  +  e2)3 


(40) 


Remark  9.  An  inspection  of  Eq.  (40)  revels  that  p*  \=  |(2  +  e2/9) 
is  a  critical  value  in  the  sense  that  #  =  tz /2  when  ft  =  /3*.  where 
we  observe  that  A(e,  ft9)  <  0  and 

Cm\p=p*  =  ~  ^  +  COS  ^  — ^  \/3  +  62  j  .  (41) 

More  interesting,  however,  is  the  following:  Since  e  1  is  as¬ 
sumed,  p*  is  very  close  to,  but  strictly  greater  than,  4/3,  the  theo¬ 
retical  value  of  p  for  monatomic  gases  under  standard  conditions; 
see,  e.g.,  Ref.  [14,  p.  80],  recalling  the  relation  p  =  (y  +  l)/2. 

Remark  10.  For  p  >  2,  it  can  be  shown  that  A(c,  p)  <  0  for 
c  <  c0,  where  e0  <  1  is  given  by 


eo  := 


27 p2  -  12 p  +44-  y/(P  -  2)  (9 P  -  10)3 


8  (P  -  1) 


(42) 


while  A(e0,P)  >  0  for  c  e  (c0,  1).  However,  in  spite  of  the  fact 
that  A(e0,  P)  =  0,  the  wave  speed  cm,  unlike  its  counterparts 
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f’*  (3  =  1.333 


Fig.  2.  /,  t fr  vs.  £  for  fj  =  4/3  and  6  =  0.1,  which  implies  cm  «  1.068,  and  Re  =  1. 
Solid:  /  vs.  £  (Eq.  (43)).  Broken:  i/  vs.  f  (Eq.  (45)). 

corresponding  to  the  LW,  Kuznetsov,  and  RSGC  equations,  never 
suffers  a  bifurcation.  This  is  easily  established  by  noting  that  A  < 
0,  for  every  e  >  0,  when  j3  e  (1,2)  and  recalling  that  Eq.  (37) 
admits  only  one  positive  root  when  >  2.  It  is  also  noteworthy, 
again  assuming  /3  e  (2,  7),  that  A(  1/2,  ft)  =  0,  where  p  \= 
5^(73  + 13^13)  «  2.220. 


Table  1 

P  =  1.047,  cm  ^  1.054,  l  «  36.319. 


Model 

0(F|n) 

d(u„) 

5(4) 

LW  (n  =  0) 

0.0798 

0.00691 

1.8857 

BLSC  (n  =  1) 

0.0798 

0.00001 

1.8857 

Kuz.  (n  =  2) 

0.0734 

0.00029 

1.6964 

RSGC  (n  =  3) 

0.1436 

0.00249 

3.7050 

Table  2 

P  =  p\Cm  «  1.068,  l  & 

28.479. 

Model 

2(F|n) 

d(u„) 

5(4) 

LW  (n  =  0) 

0.0725 

0.01265 

1.5039 

BLSC  (n  =  1) 

0.0725 

0.00059 

1.5039 

Kuz.  (n  =  2) 

0.0376 

0.00324 

0.4306 

RSGC  (n  =  3) 

0.0931 

0.00064 

2.0351 

Table  3 

P  =  2.350,  cm  «  1.116,  l  z 

a  17.000. 

Model 

2(F|n) 

d(vn) 

5(4) 

LW  (n  =  0) 

0.0297 

0.05987 

0.0215 

BLSC  (n  =  1) 

0.0297 

0.00881 

0.0215 

Kuz.  (n  =  2) 

0.1556 

0.03312 

2.6180 

RSGC  (n  =  3) 

0.0978 

0.02295 

1.6257 

Remark  11.  Since  c  1  implies  e  <  1/2,  it  follows  that  cm  e 
(1,  2)  for  every  ft  e  (1,  2)  U  (2,  7);  i.e.,  like  those  of  the  four 
weakly-nonlinear  models,  the  travelling  wave  profiles  given  below 
in  Eqs.  (43)  and  (45)  propagate  at  supersonic  speed. 


4.3.  Exact  TWSsfor  the  velocity  and  density  fields 


Returning  now  to  Eq.  (35),  separating  variables,  and  then 
integrating  using  partial  fractions,  we  obtain,  after  enforcing  the 
wavefront  condition/(0)  =  1/2,  the  following  exact,  but  implicit, 
TWS  for  the  velocity  field  (see  also  Ref.  [11,  Eq.  (3.41)]): 


5  (Cm)  : 


x  tanh 


-  In 


1  \P  —  2  —  2  eg,  QB  -  1) 
Re(cm  —  Cm)  |  VQ(Cm) 

/3-2  +  2cg,(l-f)1 

VQfcn)  J 

fcm 

,\/5 f2  -  K2C2  +  p  -  2)  +  C2  -  1 


///m 


/  6  (0,  1). 


(43) 


Here,  we  have  set/m  :=  cm/c  for  convenience  and  we  note  that 
the  shock  thickness,  /,  of  the  kinks  (i.e.,  integral  curves  of  Eq.  (35)) 
described  by  Eq.  (43)  is  given  by 


_ 2c (1  -  2/m)2 _ 

Re  1 2/3  -  ecm(l  -  2 fm)2  +  4(fm  -  1)|  ' 


(44) 


In  turn,  using  Eqs.  (33)  and  (43)  it  is  easily  shown  that  t/g  where 
we  now  let  ^(§)  :=  1  +g(§)  denote  the  dimensionless  density,  is 
given  in  terms  off  by 


VK?)  =  i  + 


Cm  -  C/(|) 


fra 

fm-m' 


(45) 


In  Fig.  2  we  have,  under  the  assumption  that  the  propagation 
medium  is  a  monatomic  gas,  plotted  /,  if  vs.  §  for  Re  =  1, 
which  we  take  without  loss  of  generality,11  and  e  =  0.1.  As  one 


11  Because  from  the  standpoint  of  numerics,  setting  Re  =  1  is  equivalent  to 
applying  the  re-scaling  £  i->-  f  /Re. 


can  clearly  see,  like  that  of  the  velocity  field,  the  density  profile 
assumes  the  form  of  a  kink,  the  shock  thickness  (ft)  of  which  is 
easily  determined  using  Eq.  (45): 


h  _  «1  -  2/m)2 

4/m(fm  -  1)  ' 

Fig.  2  also  illustrates  the  fact  that 


(46) 


0</<l<^</m(/m-l)_1  (V£eR),  (47) 

where,  in  Fig.  2,  fm  «  10.683  and /m(fm  -  l)"1  «  1.103. 


Remark  12.  Since  /m  >  2,  as  can  be  shown  using  Eq.  (47)  and 
Fig.  2,  it  follows  from  Eq.  (46)  that  h  >  l;  in  other  words,  viscous 
dissipation  evidently  has  a  greater  effect  on  the  density  than  it  does 
on  the  velocity  field. 


5.  Numerically  evaluated  metrics 

In  this  section  we  compare/contrast  the  four  weakly  nonlinear 
solutions  with/to  the  exact  result  given  in  Eq.  (43)  based  on  the 
following  three  metrics: 


d(u„)  : —  |cm  vn\ ,  and  t)(En)  : —  |/  tn\. 

Here,  D(F|n)  denotes  the  distance  w.r.t  the  L 2  norm  between  the 
kink  solution  profile  given  in  Eq.  (43)  and  the  Taylor  shock  profile 
of  the  nth  weakly-nonlinear  model,  where  Red  =  Re  is  to  be 
henceforth  understood. 

In  Tables  1-3,  we  have  computed  these  three  quantities  for  four 
different  values  of  /3.  In  every  table,  as  well  as  in  Figs.  3  and  4 
which  appear  in  the  next  section,  we  have  taken  c  =  0.1  and 
Re  =  1.  In  Table  1,  the  fluid  considered  is  Butane  (C4Hi0),  in  its 
gas/vapor  state,  at  20  °C  and  1  atm  [25];  next,  in  Table  2,  ft  = 
P9(r*  1.334)  approximately  corresponds  to  a  monatomic  gas;  and 
lastly,  in  Table  3,  f$  =  2.350  corresponds  to  liquid  sodium  at  110  °C 
and  1  atm  [13]. 

Considering  gases  and  liquids  separately,  we  seek  to  numeri¬ 
cally  identify  the  weakly-nonlinear  model,  in  the  case  of  each  fluid 
type,  that  yields  the  smallest  value  of  two  or  more  of  the  metrics. 
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f' F  J3  =  1.200 


Fig.  3.  /,  Fv s.|  for/?  =  1.200  and  e  =  0.1,  which  implies  cm  ^  1.062,  and  Re  =  1. 
(Color  available  online  only.)  Black:  Navier-Stokes  (Eq.  (43)).  Broken:  LW/BLSC 
(Eq.  (15)).  Red:  Kuznetsov  (Eq.  (21)  with  n  =  2).  Blue:  RSGC  (Eq.  (21)  with  n  =  3). 


f'F  J3  —  3. 100 


Fig.  4.  /,Fv  s.  |  for/3  =  3. 100  and  e  =  0.1,  which  implies  cm  «  1.167,  and  Re  =  1 
(Color  available  online  only.).  Black:  Navier-Stokes  (Eq.  (43)).  Broken:  LW/BLSC 
(Eq.  (15)).  Red:  Kuznetsov  (Eq.  (21)  with  n  =  2).  Blue:  RSGC  (Eq.  (21)  with  n  =  3). 

6.  Closure 

6.1.  Conclusions 

We  have,  in  the  travelling  wave  context,  compared/contrasted 
four  bi-directional,  weakly-nonlinear  acoustic  models  with/to  the 
ID  compressible  Navier-Stokes  equations  under  the  assumption 
of  homentropic  flow.  We  have  presented  exact  TWSs  for  both  the 
former  and  latter,  noted  critical  values  of  the  parameters,  and 
examined  special  cases.  Based  on  an  analysis  of  these  findings, 
in  particular,  the  numerical  results  that  have  been  presented,  we 
report  the  following: 

(i)  For  p  6  (1,  /?*),  i.e.,  in  the  case  of  gases,  there  is  not  a 
unanimous  “best  approximation”  based  on  our  metrics.  To  be 
more  specific,  while  Kuznetsov’s  equation  (n  =  2)  yields  the 
smallest  values  of  D(F  |  n)  and  d(ln),  the  BLSC  equation  (n  =  1) 
gives  the  smallest  value  of  d(un);  see  Fig.  3,  wherein  the  value 
of  /3  taken  corresponds  to  a  diatomic  gas  (e.g.,  air  at  20  °C  and 
1  atm)  [13],  along  with  Tables  1  and  2. 

(ii)  For  /3  >2  (i.e.,  for  liquids),  the  Burgers’  TWS  was  found  to 
be  equal  to  (in  terms  of  V  and  d)  or  better  than  (in  terms  of  d) 
that  of  the  BLSC  equation,  the  best  of  the  bi-directional  models 
in  the  case  of  liquids;  see  Fig.  4,  wherein  the  value  of  ^  taken 
corresponds  to  distilled  water  at  0  °C  [13],  as  well  as  Tables  3 
and  A.l. 

(iii)  As  compared  with  the  other  three  bi-directional  models,  the 
n  =  2  (resp.  n  =  1)  case  yields  the  smallest  value  of  |d£/dt| 


for  every  /3  in  the  range  corresponding  to  gases  (resp.  liquids); 
see  Fig.  1. 

(iv)  The  LW  equation  had  the  poorest  performance,  overall,  in 
terms  of  wave  speed,  with  the  value  of  d  (v0)  always  exceeding 
not  only  those  of  the  n  =  1,2,3  cases  but  also  that  of  Burgers’ 
equation  (see  Table  A.l). 

6.2.  Discussion 

Of  our  findings,  (ii)  is  clearly  the  most  surprising,  given 
that  Burgers’  equation  (see  Appendix  A)  is  actually  an  evolution 
equation  based  on  the  assumption  of  plane  unidirectional  flow.  A 
key  aspect  of  its  success,  as  shown  in  Table  A.1,  is  the  fact  that 
d(uBeq)  <  d(iq)  when  >  2,  where  it  is  noteworthy  that  iq  = 
^Beq  +  0(e2)  >  ^Beq- Of  course,  if  the  medium  in  question  is  a  liquid 
and  bi-directional  propagation  is  a  possibility  (as  in,  e.g.,  acoustic 
reflection  problems),  then  one  would  likely  be  compelled  to  set 
aside  Burgers’  equation  in  favor  of  the  BLSC  equation,  which  in  the 
context  of  TWSs  is  the  simplest  of  the  bi-directional  models. 

It  was  also  something  of  a  surprise,  given  the  results  of 
Refs.  [11,16],  and  the  fact  that  the  derivation  of  the  BLSC  model 
involves  fewer  expansions  than  that  of  Kuznetsov’s  equation,  to 
find  that  the  latter  outperformed  the  former,  w.r.t  two  of  our 
three  metrics,  when  only  gases  were  considered;  see  (i).  Flowever, 
one  finds  very  little  difference  in  both  the  solution  profiles  (see 
Table  1)  and  energy-rate  curves  (see  Fig.  1(a))  of  these  two  PDEs 
over  1  <  ^  <  1.1,  a  range  that  includes  the  /3 -value  of  not  only 
Butane  but  other  important  gases  as  well;  see,  e.g.,  Refs.  [  14,  p.  640] 
and  [26,  p.  788],  once  again  recalling  the  relation  /3  =  (y  +  l)/2. 

What  is  more,  the  fact  that  the  same  bi-directional  models  that 
performed  best  in  terms  of  the  metrics  V  and  $  also  yielded  the 
smallest  values  of  |  d£ /d t  |  (see  (iii))  suggests  that  the  minimization 
of  |d£/dt|,  over  the  ^-intervals  corresponding  to  gases  and 
liquids,  might  be  the  only  criteria  needed  to  identify/derive,  from 
among  such  models,  those  that  best  approximate  /  and  l ,  but  not 
necessarily  cm. 

And  while  the  LW  equation  exhibited  a  rather  poor  performance 
here  (see  (iv)),  as  did  its  lossless  version  in  tests  performed  by 
Christov  et  al.  [16],  it  should  be  noted  that  this  PDE  remains  a  topic 
of  interest  to  researchers  in  a  number  of  acoustics-related  fields; 
see,  e.g.,  Refs.  [27-29]  and  those  therein. 

6.3.  Possible  future  work 

We  now  mention  just  a  few  of  the  ways  in  which  the  present 
study  might  be  extended.  Perhaps  the  most  obvious  follow-on 
investigation  would  be  to  carry  out  a  travelling  wave  analysis 
using  the  un-approximated  (i.e.,  nonlinear)  ID  energy  equation 
[30,  Eq.  (14.5.16)]  and  the  non-isentropic  equation  of  state  for  a 
perfect  gas  [5,  Eq.  (5)]  in  place  of  Eqs.  (8)  and  (9),  respectively;  see 
also  Ref.  [31].  Another  possibility  is  to  take  the  constitutive  relation 
for  the  shear  stress  as  the  Ostwald-de  Waele  power-law  [24,  p.  25  ], 
a  formulation  that  would  likely  allow  for  TWSs  in  the  form  of 
compact  kinks  [32].  And  lastly,  following  Straughan  [33],  one 
could  assume  that  the  flow  of  heat  in  the  fluid  is  described 
not  by  Fourier’s  law  but  instead  by  Christov’s  [34]  refinement 
of  the  Maxwell-Cattaneo  law,  a  result  of  which  would  be  the 
introduction  of  the  fluid’s  thermal  relaxation  time  as  a  new 
parameter;  see  also  Ref.  [35]  and  those  therein. 
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Table  A.1 

For  Eq.  (A.2)  with  e  =  0.1  and  Red  =  Re  =  1. 


p 

^Beq 

®(5I0 

d(Ueeq) 

^(-^Beq) 

1.047 

1.052 

0.0798 

0.00136 

1.8857 

1.200 

1.060 

0.0773 

0.00158 

1.7178 

F 

1.067 

0.0725 

0.00163 

1.5039 

2.350 

1.118 

0.0297 

0.00193 

0.0215 

3.100 

1.155 

0.0477 

0.00814 

0.6203 
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Appendix.  Burgers’  equation 

As  shown  by  Crighton  [7],  under  the  assumption  of  unidirec¬ 
tional  flow,  Burgers’  equation  can  be  derived  from  the  BLSC  equa¬ 
tion  with  little  difficulty.12  In  the  case  of  right-running  acoustic 
waves,  and  retaining  x  and  t  as  the  independent  variables,  this  cel¬ 
ebrated  PDE  assumes  the  form  [2,7] 

2ut  +  2(1  +  ef)u)ux  =  (Red)-1!!**.  (A.l) 

Seeking  travelling  wave  solutions,  we  set  u(x,  t)  =  #(§(uBeq)) 
in  Eq.  (A.1)  and  then  impose  the  asymptotic  conditions  u  ^  1,0 
as  §  =poo.  On  solving  the  resulting  Bernoulli  equation  subject 
to  the  wavefront  condition  £(0)  =  1  /2,  we  obtain  the  well  known 
result 

5(f)  =  \  {l  -  tanh[2£0Beq)/4eq]} ,  (A.2) 

where  the  wave  speed  and  the  shock  thickness  are  respectively 
given  by  uBeq  :=  1  +  c^/2  and  £Beq  =  i\.  Here,  we  note  that 
v\  >  uBeq,  a  fact  easily  established  with  the  aid  of  Eq.  (18),  and  that, 
apart  from  the  difference  in  the  wave  speeds,  £  =  F  for  n  =  0,  1; 
i.e.,  the  plot  of  £  vs.  §  is  identical  to  the  F  vs.  $  solution  curves 
corresponding  to  Eqs.  (11)  and  (12). 

Below  in  Table  A.l,  for  the  same  values  of  j3  used  to  generate 
Tables  1-3  and  Figs.  3  and  4,  the  Burgers’  TWS  is  compared  with 
its  exact  counterpart  (Eq.  (43)),  using  once  again  the  three  metrics 
defined  in  Eqs.  (48). 
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